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ABSTRACT 

Two dimensional distortions of a panel structure due to 
welding are analyzed by the finite element displacement method 
with an assiamption of elastic deformation during welding 
process. Construction of stiffness matrices relevant to the 
welding problem is demonstrated for the one and the two 
dimensional cases. Computer programs are presented. 

Computations have been carried out using the one 
dimensional experimental value of unconstrained angular change 
along the welded edge and its equivalent constrained welding 
moment as an input to the computer program. 

These computed results are compared to the particular 
experimental values. . .Tre results indicate that an analytical 
approach using fir ement method can predict the distor- 

tion phenomena of a complicated two dimensional structure with 
a reasonable degree of accuracy. 
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Several recoiranendations are made concerning further 
investigations aimed at evaluating relevant value of angular 
change and equivalent constrained welding moment. 
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NOMENCLATURE 


a, b 


Length of a finite element in x- and y- 
direction, respectively 


C 


Equivalent welding moment for one 
dimensional deflection 


• 


Equivalent welding moment along y-direction 
as in Figure 6 


c 

Y 


Equivalent welding moment along x-direction 
as in Figure 6 


C 


Coefficient matrix 


S' 


Matrix defined in equation (49) 


D 


Matrix defined as D'/3 


E 


Young's modulus of elasticity 


F^,FT,GT,HT,qT 
— 'X''*"y ~ — .2- 


Matrix defined in Appendix A 

Transpose of the matrix F , F , G, H, and q^, 
respectively. ^ 


I 


t^/12, in.^/in. 


k 


Stiffness matrix for one rectangular finite 
element 


k , k 
—wx ~wy 


Additional stiffness matrix due to welding 
defined in equation (58b) and (58c) , 
respectively 


L , L 
X ' y 


Length of one panel structure in x- and 
y-directions , respectively 


M, N 


Number of finite element in x- and y- 
directions, respectively 


P,P ,P ,P„ ,P„ 

~ '—X '-~y ~5,x '—ily 


Matrix defined in equation (64), (65), (66) 

and (67) , respectively 


Q / Q 
~x ~y 


Load matrix due to welding defined in 
equations (58d) and (58e) , respectively 


t 


Thickness of plate 


T 


Matrix defined in Appendix A 
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1 — 1 
1 


Inverse matrix of ^ 


Ut 


Total strain energy 


u 

p 


Elastic strain energy 




Welding energy 


w 


Deflection in plate 


w 


Weight of electrode consumed, gram/ cm. 




Weight for numerical integration 


X, Y, z 


Ordinally coordinates for a finite element 


X. , Y. 
1 3 


Numerical integration points in normalized 
coordinates 


0 


Constant defined as Et^/12 (1-v^) 


V 


Poisson’s ratio 




Free joint angular change 


e . 


Constrained angular change 


6 


Variational notation 
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I INTRODUCTION 



A. Background 

Welding is used extensively in the fabrication of many 
structures, including ships, airplanes, buildings, pressure 
vessels, etc., providing many advantages over the techniques 
such as riveting, casting, and forging. However, distortion 
problems are always encountered during welding. Distortion 
often prevents the "aciiievement ..of design tolerances, reduces 
joint strength by mismatching, and imparts initial deflection 
into structural members. 

A typical structural member is a panel structure, v/hich 
is composed of a flat plate and longitudinal and transverse 
stif.feners fillet welded to the plate, as shown in Figure 1. 
Angular changes produced along the fillet welds cause out-of- 
plane distortion of the panel. The excessive out-of-plane 
distortion reduces the buckling strength of the panel that 

( 3 ) 

is subjected to compressive loading. Corrugation failures 

of bottom plating in some welded cargo ships are believed to 
be due primarily to the reduction of buckling strength of the 
plating with excessive initial distortion. 

Welding distortion is very complicated phenomena 
associated with many fac-tors;. 

1. Heat conduci:.iv I !-y of weldment which controls the 
distribution of temperature gradient of welded metal. 

2. Thermal expansion coefficient which controls the 
expansion and shrinkage of welded metal. 
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Figure 1 Panel Structure with Longitudinal and Transverse Stiffeners 
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3. Moduli of elasticity which governs the rigidity 
against deformation. 

4. Yielding strength which governs the size of the 
plastic zone near the welds. 

5. Degree of constraint which governs the freedom of 
movement during weld deformation. 

6. Amount of. heat input used which governs the size of 
the melted zone. 

7. Amount of weldment used. 

Therefore, it is extremely difficult to analyze the 

welding deformation phenomena by a pure analytical way. Yet, 

only one dimensional deformation of panel structures has been 

analyzed with the assumptions of elastic deformation during 

( 1 ) 

welding process. Obviously, welding deformations are 

neither pure elastic nor pure plastic phenomena, but the 
plastic zone near the weldments is so small compared to panel 
structural size that it is believed to be closer to the 
elastic behavior. 

Therefore, the analytical approach of the one dimensional 
deformation prediction can be treated as an elsstic deforma- 
tion by considering moduli of elasticity, degree of constraint, 
and amount of weldment used. 

B. Previous Investigations 

(1 3) 

Masubuchi, et. al. ' have investigated the one 
dimensional plate deformation during welding with a minimxim 
strain energy concept using elastic bending theory. 
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a. Free Joint 



llKjld 3E/iM 




Figure 2 Distortic/ii; Caused by Angular Change in Two Types 
of Fillet tVelded Structures 
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If a fillet joint is free from external constraint, the 
joint simply bends to a polynomical form having a knuckle 
at the weld as shown in Figure 2a. However, if the joint is 
constrained by some means, a different type of distortion is 
produced. For example, if the stiffeners are welded to a 
rigid beam, as shown in Figure 2b, the angular changes at 
fillet welds cause wavy or arc-form distortion. 

In the simplest case in which sizes of all welds are the 
same, the distortions of all spans are equal and the distor- 
tion, w, which can be expressed as follows; 

6/Jl = [1/4 - (x/t - 1 / 2 ) 2 ] .0 ( 1 ) 

where 6 = deflection of plate 

9^ = angular change at a fillet weld with constraint 
' I = length of span 



The amount of angular change, 9., in constrained 

structures is smaller than that in a free joint, 9^. This 

indicates that a certain amount of energy is necessary to 

decrease the angular change from 9^ to 9. If the necessary 

energy is represented by U , it may be given by the following 

w 

equation. 



9o-e 



dU 



w 



d(9^-9) 



(2) 



On the other hand, the strain energy stored in the 
constrained plate, U^, can be expressed using the elastic 
beam theory . 



U = 



E'l 



Z 



(3) 
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where E' = E/1 - v‘ 



I = tVl2 



Since Up increases but on the contrary decreases as 

the constrained angle e increases, the condition for equili- 

brixam of this system requires that the total strain energy 

U = U + U should be minimum. Furthermore, for the sim~ 

"C w p 

plicity of the problem, the ratio of incremental welding 
energy change to angular change is assumed to be linear as in 
the following equation. 



dU 



w 



aTe^) = 



(4) 



where C = a coefficient determined by the weight of the 
deposited metal and by the welding procedures. 
From equations (2) and (4) , welding energy per unit 
length of width can be expressed in terms of C, 9, and 9^, 



»w = I 



(5) 



Accordingly, the condition of equilibrium is as follows: 



3U , 

c(6^-e) + 0 



( 6 ) 



From equation (6) , relation betv/een 9 and 9^ can be 



expressed , 



9 = 



1 + ^/C 



(7) 



The values of C can be determined by using relationship 
obtained by experimentally measuring 9^ and 9 . Empirical 
formula for the calculation of C has been proposed by 



17 



Thickness of plate, in. 




Figure 3 Variation of Angular Change of Free Fillet Welds, 
0Q, As a Function of Plate Thickness, t, and 
Weight of Electrode Consumed per Weld Length, W. 
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. (3) 

Masubuchi. 

C = t"^/(l+W/5) kg-mm/mm radian (8) 

where W = weight of electrode consumed, gram/ cm. 
t = thickness of plate, mm. 

The experimental results of determining the unconstrained 

angular change as a function of plate thickness, t (mm) , and 

weight of electrode consximed per weld length, W(gr/cm), have 

( 2 ) 

been given by Hirai , et. al. as in Figure 3. 

Therefore, for a given plate geometry and welding con 
condition, the deflection of the one dimensional case can 
easily be estimated by using equations (1) , (7) , (8) , 

and Figure 3. 

C. Aim and Purpose of Present Studies 

As stated previously, welding deformation has been 
analyzed only for the one dimensional case, but in reality 
most of the structural members are encountered in two 
dimensional deformation, in which the deformation, w, varies 
along the x- and y-directions as well as the unconstrained 
angular change, 0, and equivalent constrained welding moment, 
C. With the idea of previous investigations which state that 
the total strain energy in the plate during welding deforma- 
tion becomes minimum, it is possible to extend the one 
dimensional deformation analysis to the two dimensional case. 
However, the simple mathematical solution using plate bending 
equation is rather difficult due to the uncertainty of the 
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equivalent load and the boundary conditions relevant to the 
welding deformation problem. 

The difficulties can be easily overcome by formulating 
the problem into a finite element model because it can be 
self-adjusted to take care of the equivalent load and equiva- 
lent sitffness matrix due to welding. Furthermore, in this 
way the relation of deformation phenomena to the constrained 
angular change, 9, and the equivalent constrained welding 
moment C can be easily seen. 

Therefore, it is intended to handle the two dimensional 
deformation with the finite element displacement technique. 
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II FORMULATION OF FINITE ELEMENT EQUATION 
A. One Dimensional Case 

In this study, computer programs are programmed only 
for the two dimensional finite element analysis, but the 
one dimensional finite element fo3rmulations are presented to 
visualize the inside of the complex matrix assembling and 
computational procedure; furthermore, to compare the results 
to that of the previous investigators. 

For the one dimensional case, as shown in Figure 2b, if 
deflections are elastic and small, the deflection can be cal- 
culated using elastic theory when the welding equivalent loads 
and the boundary conditions are prescribed. Even though the 
deformations are assxamed to be elastic for the welded struc- 
tural case, the loading and boundary conditions to be applied 
for this problem are not clear. 

As introduced by Masubuchi , et. al. , if the concept 
of equivalent constrained welding moment, C(lb-in/in rad.) 
is being used, the uncertainty of loading conditions may be 
overcome; for the boundary condition it may be safe to assxame 
to be simply supported for which only the angular change along 
the welded edge can be allowed. Furthermore, for simplicity 
of the problem, the deflections are considered only for the 
case of one span length of a simply supported beam to elimi- 
nate the complexity of interaction due to the statically 
indeterminate structural effects, as in Figure 4a. 




if 




b. Local coordinate and nodal points for a 
finite beam element 



Figure 4 One-span Beam Element and the Finite 
Element Local Coordinate 



22 



Fox the one dimensional plane stress, the strain energy 
can be expressed simply by: 



U = 1/2 Iff e a dxdydz 

p / •> ■> XX ^ 



(9) 



Using the stress-strain and strain-curvature relation for the 
beam element, 



a = Ee 

X X 



e = -z 

X 



d^w 

^ 2 
dx 



(10a) 

(10b) 



Therefore, the strain energy can be expressed in terms of 

deflection and rigidity of the plate: 

2 

dx (11) 



p 2 



^2 
d w 



where 



0 

t/2 b 



dx"” 



1 = 2 / 
0 



- 2, , bt' 

/ z dzdy = yt" 

0 



On the other hand, the welding energy can be expressed 

in terms of 0^, 9, and C, as in equation (5). For this 

problem it is more convinient to express the welding energy 

in terms of deflection rather than the angular changes. From 

equation (5) , 

^ ^ c 2 

^ 0 0 ^° 

the constrained angular change can be expressed by the 
equation 

e = aw/ax (12) 

Therefore, the welding energy expressed, 



( 13 ) 
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T 7 be r / r\ ^ \ -3 

W ^ Q O o X 



where the values of C and 0 are constant which will be 

o 

determined by the particular welding condition and the 
geometry . 

The equilibrium condition for the system requires that 
the total strain energy has to be minimized, in other words, 
the variation of the total strain energy has to be zero, 
which becomes 

6U„ = 6U^ + 6U , = 0 (14) 

i p w . 

using equations (11) and (13) , equation (14) becomes 



El ^ 

2 i 


f,2 ' 

d w 
^ 2 


2 

dx 


+ 6 


be f 
2 n 


dw. 2 , 

(0 - T— ) dx 
0 dx' 


_ 0 


[dx 






_ 0 





But the variation of the constant term is zero, therefore. 



a 

El / 
0 






be ^ 

2 ^ 
0 



,d6wv ,dw. , 

(“^H~) (j— ) dx - 
dx dx 



a 

be / 
0 



,d6w, , 



(15) 



For a beam element as shown in Figure 4b, the deflection 

can be assumed to be a general cubic function of x, 

2 3 

w = Cj^ + C 2 X + c^x + c^x (16) 

and when the nodal displacements and the nodal angles are 
the unknown variables which are to be calcualted, the coeffi- 
cients (Cj^, C 2 , c^/ c^) in equation (16) can exclusively be 
expressed in terms of the nodal displacements and angles 
considering the requirements of continuity at the nodal points 
w(x=0) = w^ = Cji^ 
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^(x=0) = 

2 3 

w(x=a) = W2 = + c^a + c^a + c^a 



•^(x=a) = ^2 '*' ^^3^ ^*^4^ 



The above requirements can be expressed in the matrix form, 



q = T C 



(17) 



where 



T = 



10 0 0 



0 10 0 

, 2 3 

1 a a a 



0 1 2a 3a^ 



(17a) 



C = 



(17b) 



and 



q = 



w. 



w 



X, 



w. 



w 



X, 



(17c) 



From equation (17) , C can be expressed 
C = t”^ q 



(18) 



If the new matrix G is defined such that 
G = t“^ 



(19) 



then, 



C = G q 



( 20 ) 




IM 
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where 



G = 



1 

0 

-3/a' 



0 

1 

•2/a 



0 

0 

3/a 



0 

0 

-1/a 



■2/a^ 1/a^ 



(21) 



2/a'' 1/a 

The second derivative of deflection function becomes, 

^ = 2c 6c X 
2 3 4 

dx 

which is in matrix form 



^2 
d w 

2 

dx 



= H C 



where H is the matrix defined by 
H = (0 0 2 6x) 

then, 

= 6 (H C) (H C) = 



( 22 ) 



(23) 



d^6w 




fn2 

d w 


^ 2 
dx 

\ J 




2 

[dx J 



('Sc'^) (A C) 



but, using the equation (20) , 



d^6w 


r^2 >1 

d w 


[dx J 


^ 2 
[dx 



= (6q^) (G^ H G q) 



(24) 



The matrix expression for the first derivative of the 

deflection w can be expressed in the same manner, 

dw , ^ _ 2 

^ = C2 + 2C3X 3 c^x 

using and C 
dw 



dx 



= F C 



(25) 
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where 

P = (0 1 2x 3x^) 

Finally, using equation (20) 

^ = («/ ) w'' 



and 



d6w' 


dw 


[dx J 


dx^ 



m mm 

= (6q ) (G^ F F G q) 



(26) 



(27) 



(28) 



From equations (24) , (27) and (28) , the minimum energy 

equation (15) can be written in the matrix form, 

^ m m m 

El / (6q^) (G^ K G q)dx 

0 ~ 

+bc / (6q'^) (g"^ f'^'^ F G q)dx 
0 " " 



-be / (6q'^) (g"^ F*^)dx = 0 (29) 

0 ^ ^ " 

. T 

Noticing that the matrix q , q , G , and is not the 

function of x, but only of the function of dimension of the 
plate, and displacement and angle at the nodal points; there- 
fore, equation (29) can be simplified to 

{El g"^ [/ (h"^ H)dx]G q 

r ^ 0 ~ ^ 

+bc g"^ [ / (f"^ £) dx] G q 

~ 0 ^ 

-be G^ [/ (F'^)dx]} = 0 
0 ^ 



( 30 ) 
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in which the only integrations necessary to be carried out 
are the bracketed quantities/ which become after integration. 







0 


0 


0 


0 


a „ 




0 


0 


0 


0 


/(H^ 


H)dx = 










0 ~ 




0 


0 


4a 


6a^ 






.0 


0 


6a2 


12a^ 






f 

\J 


V 


0 


0 






0 


1 


2a 


3a 


/(F^ 


F) dx = 










0 




0 


2a 


4a 


^ 6a 






.0 


3a 


2 ^ 
6a 


^ 9a 



^ T 
/(F^)dx 

0 



'o 

1 

2a 

^3a 



2 



(31) 



(32) 



(33) 



Equation (30) is the final equation for one element 
formulated by the finite element displacement method. The 
first quantity is defined as a stiffness matrix for an 
element, the second is the additional stiffness matrix due 
to welding, and the third is the equivalent load due to weld- 
ing. They are as follows: 
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Stiffness matrix for an element, 



= El 



T ^ T 

G [/CH 
0 

_ 


H) dx j G 




'' 


12/a"' 


6/a^ 


-12/a^ 


6/a^ 


6/a^ 


4/a 


-6/a^ 


2/a 


-12/a^ 


-6/a^ 


12/a^ 


-6/a^ 


6/a^ 


2/a 


-6/a^ 


4/a ^ 



(34) 



Additional stiffness matrix due to welding, 
k„ = be F)dx]G 

'O 0 0 0 ' 

0 0 0 0 

0 0 0 0 

.0 0 0 bcj 



(35) 



Load matrix due to welding , 

Q = be G"*^ [/ F"^ dx] 

^ 0 ~ 

0 

= ° (36) 

0 

be© 

i oJ 

At this stage, assembling of equation (30) has to be 
eonsidered when the one-span beam is divided by many finite 
elements, but the one-span beam element of Figure 4a is 
ehosen to be two finite elements as in Figure 4b; then, only 
element 1 can be used to predict the deflection and angular 
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change at the nodal points by symmetric argument, 
case, the simultaneous equations to be solved are 
equation (30)), 

'I + “x, - I ”2 + “x = ° 

a 1 2 



In this 
(from 



(| - I “2 * ”x2>’ ' “ 



r 6 El / 2 ,2 \ 1 A 

Swj <- S «1 - «x, a «2 - “x >' = “ 

^ 12 



<1 ''1 + - I «2 <2 + i|)«x2 - = “ 

(37a-d) 

From the simply supported boundary condition at node 2 , 



W 2 =' 0 ; therefore, 6 W 2 = 0 . 

From symmetry condition at node 1, w =0; therefore, 

^1 

6 w = 0 . 

^1 



Furthermore, the displacement at node 1 and the angle at 
node 2 is not known, which means that the variations of these 
variables are not zero; therefore, the only way to satisfy 
equation (37a) and (37d) is the bracketed value in each equa- 
tion has to be zero, which leads to final equation for w 

and w , 

^2 



— w, + w =0 
a 1 Xj 



I „1 + (2 * 



0 



(38a-b) 
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From equation (38a-b) , displacement at the middle point and 
angular change at the edge expressed as, 

-a6o 



^1 2(1 + El/abc) 



6 ^ 



^x^ 2 (1 + El/abc) 



(39) 



(40) 



It is worth to compare the results of equations (6) and 
(7) with equations (39) and (40) . In equations (39) and (40) , 

El _ btVl2 _ t^/12 



abc 



abc 



ac 



and a = Jl/2 is to be used. Then, 



w, = 



-ae 



n 



1 4(1 + 2EI/£c) 

_ ®o 

*'^X 2 “ 1 + 2El/ilc 



(41) 

(42) 



f 1 6 ) 

where I is defined by unit width basis which is 
I = t^/12. 

Equation (42) agrees with equation (7) , as also does 
equation (41), if x = 1/2 and 6 of equation (7) is plugged 
into equation (1) . The negative sign of equation (41) means 
that the coordinate and the positive angle at the edge has to 
be chosen, as in Figure 4a. 

As shown, the results from the finite element method is 

(1 3) 

exactly the same expression proposed by Masubuchi, et.al. ' 
for the one dimensional deformation, and in this way, the 
equivalent loads and the stiffness matrix due to welding 
become easily visualized. 



m 
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B. Two Dimensional Case 

1 . Formulation of the Problem 

In the two dimensional analysis, all procedures are the 

same as that of the one dimensional analysis except that the 

integrations are carried out by the numerical methods and the 

methods of assembling are mentioned. The deflection may be 

(12 14) 

approximated by the bi-cubic function for one element, ' 
as in Figure 5. 

2 2 3 

w = c^ + C 2 X + c^y + c^xy + c^x + Cgy + c^x 



2 2 3 3 2 2 

+ CgX y + CgXy + c^^^y + c^^^x y + c^2^ ^ 




(43a) 






(43b) 





•^•5 22^- 32 

+ 3c^^x y + 3c^gX y 



(43c) 





3x3y 




(43d) 



32 



i 



13 = 






‘U = UJ 3 



^(4 =. (A/y4 
^IB - l^y4 
^ ‘Lie-uTfc^^ 



%iO = 1^X3 

^ 11 = 3 

^IZ= UJx^3 




= U7| 
^2/= UJxi 

= Ujg , 

^4 = I 



^6 = ^-XX 

ig= XAT;t^2. 



Figure 5 Local Coordinate of a Two dimensional Rectangular 

Finite Element and Generalized Displacement Numbering 
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From equation (43a-d) , the matrix equation can be defined 
such that at each nodal point, the generalized displacements 
(deflections, angles in x and y, and twisting angles) are 
expressed in terms of the geometry of the element and the some 
unknown coefficients which is the form of, 

q = T C (44) 



where ^ is 16 by 16 matrix, £ is 16 by 1 matrix, and q is 
16 by 1 matrix defined as in Appendix A. 

From equation (44) , the unknown coefficient matrix, C, 
can be expressed in terms of geometry and the generalized 
displacements , 



C = G ^ (45) 

where 

. G = t"^ (46) 



However, the actual calculation of the inverse of T is 
carried out by the computer in this analysis. 

From the elastic plate theory , ^ the strain energy 

stored in the plate during the elastic deformation is ex- 
pressed by _ 



a b 

= § / / 

P ^00 



2 2 

3 w , 3^w 

2 ^ 2 
3x 3y_ 



2(1-V) 



.2 .2 

3 w. 9 w 

2 2 
3x 9y 



3^w 

3x3y 



2 



dxdy 



This can be written in the matrix form, 
a b 



( 7 ) 



which is 



U = 
P 



1 

2 



0 0 



D K dydx 



(47) 



where 
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K = - 



9 w 
2 

9x 

^2 
9 w 

c. 2 

9y 

,9 2w 



'9x9y 



= -H C 



(48) 



in which H is given by Appendix A, 



D' = ^ 



12(l-v^) 



in which 



B = 



Et 

12(l-v^) 



Cl V 0 

! 

V 1 0 

1-v 



0 0 



= 3 



and ' 



fl V 0 1 



D = 



V 

0 



1 

0 



0 

1-v 

2 



(49) 



Therefore, equation (48) becomes after taking the 
variation with respect to generalized displacement, 



6U = / / (6K )B D K dydx 

P 0 0 



= / / (6C ) (H 8 D H Odydx 
0 0 



T a b pp 

= 6q 3 ^ [/ / IP D H dydx]G q 
0 0 ~ ~ 



( 50 ) 
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The welding energy for the elements shown in Figure 6b 
can be written as, 



Therefore , 



C„ a 



^w “ ^ ^®o ■ ^ 






3Wi j 



(51) 



which becomes , 



""w = 



/ (^) dx + -; 5 - / (— ) dy 



2 0 



C 9 / ^ dx + C 0 / dy 

yo^ay xoq9x^ 



i(C + C )6^ 
2 y X o 



(52) 



and taking the variation, 

a 



= 

w 



c . / 



9 6w 



9yJ 



3w 

Wyj 



b r 



dx + C / 
^ 0 



K> 

CO 


faw' 




axj 


19yJ 


dy 



C 0 / 

^ ° 0 



a^iw 



ayj 



dx + C 0 / fe- 

X o Q (ax 



dy 



(53) 



in which the first bracketed term is the additional stiffness 
matrix due to welding, and the second is equivalent load 
matrix for the welding. In matrix form, equation (53) can be 
written, if some matrices defined as below; 




C 




aw 

ax 



(54a) 
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3w 



F 






C 



C 



T 




(54b) 



where the subscripts x and y denote that the integration 
should be carried out along the x- and y-direction, respec- 
tively, and the matrix F^ and F^ are defined in Appendix A. 
Therefore , 



9 6w 


'9w' 


9yJ 


[syj 


96w 


'9wl 


9x^ 


,9xJ 



( 8 ^) (F^ F^ C) = Fy F^ G (55a) 



(Sc'^) (f"^ F C) = (6g'^)G'^ f"^ F G g (55b) 

^ ~X '-X ^ ^ ^X '-x ^ 



Using equations (54a-b) and (55a-b) , the welding term 
can be written as: 



6U^ = 6q’^ g"^ [/ (F^ F )dx]G q 

+ 6q'^ g'^ [/ (F^ F )dy]G q 

- 6q'^ C 0 g’^[/ (F^)dxl 

~ y o ~ Q -X 

- 6q'^ C 0 G^[/ (F^)dyJ 



(56) 



As in the one dimensional case already shown, the total 
variation should be zero to be in equilibrium of this system. 
Therefore, 



6U 



T 



6U + 6U = 0 
p w 



which is from equations (50) and (56) , 
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+ C 0 / dyjl = 0 (57) 

X o Q -y ^ j 



These are the 16 by 16 simultaneous equations which have 
to be solved if only one element is concerned as shown in the 
one dimensional case, with the proper boundary and symmetry 
conditions. From equation (57) , the stiffness matrix for a 
plate element, 

k=G^3//H^DH dydx G (58a) 



The additional stiffness matrix due to welding along the 
x-direction of the plate. 



T ^ T 

k = G C / F F dx G 
-wx ~ y ^ --X ~x ~ 



(58b) 



The additional stiffness matrix due to welding along the 
y-direction of the plate. 





C 



X 





dy G 



(58c) 



The equivalent load matrix due to welding along the 
x-direction of the plate, 

Q . = g"^ C 0 / f'^'^ dx (58d) 

- y o 0 



The equivalent load matrix due to welding along the 
y-direction of the plate. 



T ^ T 

Q = G C 0 / F dy 

Xy -- X o Q -y 



(58e) 
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Figure 6 Various Finite Element Having Welded Edge 
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2 . Assembling Procedure 

For all the finite elements except those having the 
welded edge, the stiffness matrix necessary to be assembled 
is the expression given by equation (58b) , and for the finite 
elements having welded edges, three cases are involved. The 
first case is elements having welded edges along the x-direction 
as in Figure 6b; the second case is along the y-direction, as 
in Figure 6c; and the thi'^d case is along x- and y-directions 
as in Figure 6b. The asseiublijig procedure is the expression 
of equation (57) for the overall plate which is to be solved, 
which is described in FEABL user's manual. Also, from 

the symmetry condition, which -will be described in the follow- 
ing section, only one-quarter part of the plate is actually 
being considered. 

Example of the assembling procedures for the two finite 
elements is given in Appendix B. 

3. Symmetry and Boundary Conditions 

After assembling, relevant symmetry and boundary 
conditions have to be applied to solve the simultaneous 
equations, as shown in the one dimensional case. Along the 
welded edge of x and y, the simply supported edge condition 
requires that all the displacements at that nodal point have 
to be zero, and for the x-directed edge, the angular changes 
in the x“<3i^ection are not allowed, In just the same way, 
for y-directed edge, the angular changes in y are not allowed. 
The symmetry condition requires that the first derivative 
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along the symmetric axis should be zero. For example, the 
symmetry and boundary conditions for the case given by 
Figure 7 are as follows ; 

Boundary conditions : 

At node 6/ 12, 18, 24, 30, 36, 42, 48, 54 
w = 0 
= 0 

At node 4 $^ bl , 52, 53, 54 

w = 0 




Symmetry conditions : 

At node 1 , 2 , 3 , 4 , 5 , 6 




At node 1, 7, 13, 19, 25, 31, 37, 43 
At node 1, 6, 49 

4 . Nvimerical Integration Method 

Usually in the finite element analysis, the integrations 
involved are carried out by the numerical methods, which can 
be easily calculated using the computer. For the present 
analysis, the nine Gaussian points are chosen for the rectan- 
gular finite element , as shown in Figure 8c. In actual 
integration, it is very convenient to use the normalized coor- 
dinate system, as will be described below. From Figure 8a and 
8b, normalized coordinate system can be defined as 
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X = |-(X + 1) 
y = |(Y + 1) 

and 

dx = ^ dX 
dy = ^ dY 

Therefore, the arbitrary integration 
a b 

f(a,b) = / / P(x,y)dxdy 
0 0 

becomes in the normalized coordinates, 

,11 , 
f(a,b) = 1^ / / P[|(x+1), f(Y+l)] 

^ -1 -1 ^ ^ 

^ Equation (62) can be integrated by numerically 
9-Gaussian points, which become 
,33 

f(a,b) = ^ E I W.W. P(x,y) 
i=l j=l ^ 

where 



X = 


l^^i 


+ 1) 






y = 


|(Yj 


+ 1 ) 






W. : 
1 


Wi 


= 5/9, W 2 


= 8/9, 


= 5/9 


W . : 
D 


Wi 


= 5/9, 


= 8/9, 


= 5/9 


X. : 
1 


^1 


= -/TVS, 


X 2 = 0 , 


, X 3 = /IV5 


^3 = 


^1 


= -/I5/5, 


0 

II 

CN 


, Y3 = /I5/5 



(59 a) 
(59b) 

(60a) 

(60b) 

(61) 

(62) 

using 

(63) 

(63a) 

(63b) 

(63c) 

(63d) 

(63e) 



(63f) 
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Figure 8 Ordinally and Normalized Coordinates for an 
Element, and 9 Gausian Points for 
Numerical Integration 
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Therefore, using the concept of normalizing and 
numerical integration described above, the integrations 
involved in the present analysis such as equation (58a-e) 
can easily be written, which are programmed in Appendix D. 
From equation (58a) 
ab fi T ^ ^ 

k= {^) G-^ [ E E W.W. P(x,y)]G (64) 

~ i=l j=i 1 J 

where 

P(x,y) = D H 

X, y, Wj , and are defined by equation (63a-f ) . 

From equation (58b) , 
aC 3 

G i ^ W. P^(x,y)]G (65) 

^wx Z — i=l ^ ^ 



where 



P (x,y) = F^ F 

^ --X 



X = |(X. + 1) 
y = b 

and are defined by equation (63c) and (63e) , respectively. 
From equation (58c) , 

’^‘=x T ^ 

>Sv/y = — £ Wj Py(^-y)l£ 

where 

P^(x,y) = F-^ F^ 

X = a 

y - |(Yj + 1) 



( 66 ) 
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W. and X. are defined by equations (63d) and (63f) , 
3 3 

respectively . 

From equation (58d ) , 

T ^ 

&C “ -r^ a i.s 

1=1 

where 






T 



^ + 1) 
y = b 



(67) 



and are defined by equation (63c) and (63 ) , 

respectively . 

From equation (58e) , 

be C 3 

Qy = Wj 

where 

= Ey 

X = a 



( 68 ) 



y = |(Y. + 1) 

Wj and Yj are defined by equations (63d) and (63f) . 

As noticed above , numerical integration using the 
normalized coordinate can simplified the computing procedure 
by factoring out a and b. 
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5 . Computer Programs 

Computations of the ordinally, additional welding, and 
equivalent welding load matrix of equation (58a-e) have been 
prograjnmed by using the numerical integration technique with 
normalized coordinates as listed in Appendix D. In addition 
to these, FEABL Programs in assembling of the relevant 

matrices and solving the simultaneous equations with proper 
boundary and symmetry conditions are also used and listed in 
Appendix D. 

The Programs listed in Appendix D are only relevant to 
the panel structure welding with rectangular finite elements. 
Further, the nodal nxambering system should follow the way 
given by the sample analysis, which means that the numbering 
should start from origin to y-direction, as shown in Figure 7. 

Input data for this analysis are the overall dimensions 
of one panel structure in the x- and y-direction, thickness 
of the plate. Young's modulus of elasticity, Poisson's ratio, 
number of elements divided in x~ and y-direction, equivalent 
welding moment in x- and y-direction, and free joint angular 
change at the edge. The input data cards are described in 
Appendix C. 

Outputs are the displacements , angular changes in x- and 
y-directions , and the twisting angles at each nodal point, as 
shown in Appendix E. 

One thing to be noticed here is that for a quarter part 
of the plate, as shown in Figure 9, the proper boundary and 
symmetry conditions at the nodal points are to be generated 
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by subroutine HOLD given in Appendix D; therefore, as long as 
boundary conditions of simply supported and symmetry conditions 
as described are concerned, it is not necessary to be consi- 
dered in the input. Therefore, no matter how the number of 
finite elements are to be chosen, the only input data will be 
as described in equation C69a-h) . Furthermore, to use this 
program, the coordinate and a quarter part of the plate should 
be chosen as the same way as shown in Figure 9 . 

The units used in the sample program are pounds, inches, 
and radians. As results, for the output data, the deflections 
are in inches, the first derivatives (angular changes) are in 
radians, and the second derivatives with respect to x and y 
(twisting angles) are also in radians. 

• Program language used is FORTRAN IV, and the computations 
are carried out by an IBM 370/155 system. 
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Figure 9 Global Coordinate System and Shape of Welded Deflection 
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Figure 10 A Panel Structure Analyzed by Experiment and 
Finite Element Method 
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III RESULTS 

As computer inputs, one steel panel structural dimension, 
as shown in Figure 10, for which the experiment conducted and 
weight of consumed welding rods are given, 3 gr/cm. , 
are used to predict the deformations , angular changes , and 
twisting angles at the nodal points. Here, the value of the 
one dimensional free joint angular change along the welded 
edge, 6^, has been read frc"* the results of Figure 3 and the 
equivalent constrained welding moment, C, is being evaluated 
using equation (8) . The mamerical values of the input data 
are as follows; 

Half-length of the plate in the x- and y-directions , 
respectively, L = 16 inches and L =10 inches. 

Number of finite elements in the x- and y-directions, 
respectively, M = 8 and N = 5. 

Equivalent constrained v/elding moment from equation (8) 

c = - (25.4 X 3/8)^ _ kg-mm 

1 + W/5 1 + 3/5 mm rad. 

where W is given by 3 gr/cm from experiment. 

Therefore , 

C = 1.13 X 10^ Ib-in/in rad. 

Poisson's ratio: 

V = 0 . 3 

Young ' s modulus ; 

E = 30 X 10^ psi. 
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Thickness of the plate: 

t = 3/8 inch. (69 a-h) 

The results of the analysis using the above data are 
listed in Appendix E, and these results are plotted in the 
Figures 11 and 12 to be compared with experimental results. 

As can be noticed from Figures 11 and 12, the finite element 
results are lower value than the experiments, and the maximxam 
deflection at the midpoint of the plate which is the most 
important deflection in reality is only a half of that from 
the experiment. 

To visualize the behavior of the deformations as 

functions of 0^ and C, series of computations are conducted 

with different combinations of 6 and C as shown in Table 1. 

o 

Yet remaining, the rest of the input data are unchanged, and 
the maximxam deflections at nodal point 1 are plotted in terms 
of 0Q and C as in Figures 13 and 14 , and the maximum angular 
changes at node 6 and 49 are plotted in the Figures 15 and 16. 
Only the maximum deflections at node 1, and the maximum 
angular changes at node 6 and 49 from the results of series 
computations are tabulated in Tables 2 and 3 . 

-3 

Finally, computed results using 0^ = 55 x 10 and 
C = 4.2 X 10^ are listed in Appendix E, and also plotted in 
Figures 11 and 12 to compare with the experiment. As can be 
seen from Figures 11 and 12, these i^esults very close to 

the experiment and the reason of choosing the above 0^ and C 
will be discussed in the following section. 



SKpBZfnaNT 0 C = iJ3x/0^, do- ^5X10'^ O C = ^.2x/0 ,Q,= S5xiO 
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Figure 11 Deflection Comparison with Experiment and Finite Element Analysis 
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Figure 12 Angular Changes along the Welded Boundary Compared with 
Experiment and Finite Element Analysis 



Maxirruam Deflection at Node 1 in 



54 



to 

<u 

o 

H 




Angular Changes (9w/3y) at Node 



55 



VO 



g 

•H 

X 




Angular Changes (9w/9x) at Node 49 
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IV DISCUSSION OF RESULTS 

As has been mentioned in the previous section, the 

. 4 -3 

results using C = 1.13 x 10 and 0^ = 45 x lo are only a 

half of the experimentally measured values as shown in 
Figures 11 and 12. The reasons of this difference are not 
yet known. These may result from the assumption of the elas- 
tic behavior of the welded deformation, or from the uncer- 
tainty of the experimental measurement itself because the 

deflections and angular changes which are only of the order 

-1 -2 

of magnitude of 10 inches and 10 radians, respectively. 
Yet, other possibilities are from the value of 0^ and C. 

To investigate the behavior of the defonriation with 
respect to 0^ and C, series of computations have been con- 
ducted. The plotting of maximvim deflections at node 1 
versus C with changing 0^ as in Figure 13 shov;s some insight 
of deflection behavior with varying C. In here, it is ob- 
vious that the deflections are more rapidly changed with 
changing 0^ than C changes, which means that the deflections 
are more strongly dependent function of 0^ rather than C. 
Furthermore, the deflection reaches rapidly to an asymptotic 
value as C increases, which can be seen clearly from the one 
dimensional deflection, equation (41) where the term 2EI/ ic 
becomes small compared to 1 as C increases; therefore, in the 
limit case, the deflection can be expressed only by 

w = Zd /4. 
o 
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This behavior of the deflection with changing C are 
clearly shown in Figure 13. Also from Figure 14, in which the 
maximum deflections are plotted with respect to 0^, it is 
obvious that the deflections are rapidly increased as 6^ 
increases. 

If assumed that the difference between experiment and 
computed results are caused from the possibility of different 
values of 0^ and C, then the question is how to choose 0^ and 
C, such that the computed maximum deflection becomes the same 
to that from the experiment. From Figures 12 and 13, the com- 
binations of 0^ and C to satisfy the given value of maximum 
deflection from the experiment can be predicted. However, 
the best combination is not known, but for the known experi- 
ment value of maximum deflection it may be safe to choose 0^ 
and C such that the changes of both 0^ and C from the cal- 
culated value be minimized. In other words, superposing 
Figures 13 and 14 as shown in Figure 17, the point which 
gives the same reading of 0^ and C in both ways as shown by 

-3 

the small circle, which is read as 0 = 55 x 10 and 

o 

4 

C = 4.2 X 10 , will be the optimal combination to satisfy 
the given maximum deflection. The comparison of the results 
using this 0^ and C to the experiments , as in Figures 11 and 
12 , shows that this approach gives reasonably good results . 

As has been seen from above, there exists many 
uncertainties for the two dimensional welding analysis, but 
the welding deformation phenomena can be predicted with a 
reasonable degree of accuracy by the finite element analysis. 
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Furthermore, the empirical formula of C given by equation (8) 
may not be simply applied to the two dimensional case, and 
also does the value of 0^. If it is so, measuring only the 
maximum deflection in future experiments, the best combination 
of 6^ and C can be predicted as discussed. 



65 



V CONCLUSIONS AND RECOMMENDATIONS 
A. Conclusions 

From the discussion of the previous section, it can be 
concluded as ; 

1. The finite element analysis for the two dimensional 
welding deformations which are assumed to be elas- 
tic deformation gives reasonable prediction of 
deflections . 

2. The searching technique of optimal combination of 

0^ and C for a given deflection can be used as shown 
already, if there exists the uncertainty of 0^ and C. 

3. In the future experiments, this analysis gives the 
order of magnitude feeling of the deflection mea- 
sured after welding. It is extremely difficult to 
measure correct deflections ' in experiments because 
the order of magnitude to be measured is very small; 
therefore, the inherent uncertainties always exist. 

4. The magnitude of importance of the welding deformation 
variables such as 0^ and C are clearly seen, there- 
fore, more care should be concentrated in measuring 

0Q rather than C. 

5. Simple empirical formula for C given by equation (8) 
may not cover the wide range of the values of C; 
furthermore, as can be derived from equation (41) it 
may be logical to express C in terms of maximxam de- 
flection, free joint angular change, and rigidity. 
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such as : 

C = 6 /4w - 2EI/Z (70) 

o max 

in which 0^ is a function of a given welding 

condition and w is the constrained maximum deflec- 
max 

tion at midpoint of the span for the same condition 
of the welding as of 0^. Further, the elastic 
rigidity and plate dimension terms are expressed by 
the term of 2EI/£. 

B. Recommendations 

In this analysis, the assumption used is the elastic 
behavior of deformations. As stated in the conclusion, it is 
still believed reasonable to use the elastic theory, but for 
more accurate prediction, local plastic deformation phenomena 
should be included in the analysis. 

If the elastic assumption is to be used, more careful 
analysis to evaluate the value of 0^ and C are recommended in 
future research. Also, a more rational approach to express 
C is recommended. 

In future research, it is recommended that series of 
exper intents be conducted to measure the maximtim plate defor- 
mation only, and express C and 0^ in terms of welding and 
plate dimension. In this way, the computer inputs are always 
available for any prediction of deformation for the two 
dimensional panel structures . 
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APPENDIX A 



Definitions of Matrices 
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APPENDIX B 

Example of Assembling 



As an example, the assembling procedure for a quarter 
part of the plate divided by two finite element shown in 
Figure 18a will be discussed. 

For the element 1, as shown in Figure 18b, the stiffness 
matrix in a local coordinate can be written from equation 



(58a, b) , 




k + k 



-w 



X 



(A-1) 



where k. is 16 by 16 symmetric matrix. And from 
~ 1 ■ 

equation (58d) , the load matrix for the element 1 is. 



Stj^ 2x 



(A- 2) 



where Q, is 16 by 1 matrix. By the same way for element 2, 
~'^1 

the stiffness and load matrix can be written using equation 
(58a, b, c) and (58d, e) , respectively. 



k. = k + k + k 
~t ~ --WX -~wy 



= Q + Q 

~y 



(A- 3) 
(A- 4) 



where k. is 16 by 16 symmetric matrix and Q. is 16 by 1 
matrix . 

At every nodal point there are four unknowns (w, w , w , 

X y 

w ), therefore, after assembling element 1 and 2, the 
xy 

equations have to be solved with proper boundary and symmetry 
conditions are 24 by 24 simultaneous equations, which is in 
the matrix form; 



I 

I 



I 



m 



j 
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6 ^ [K q - Q] = 0 



(A- 5) 



where the matrix K is 24 by 24 symmetric matrix, Q is 24 by 1, 
. T 

q IS 24 by 1, and q is 1 by 24 matrix. 

All the quantities in equation (A-5) are defined by the 
global coordinate, therefore, the assembling procedure is to 
locate the quantities defined by equation (A-1, A-2, A-3, A-4) 
to the proper position in equation (A-5) . This operation can 
easily be established comparing each nodal points numbered in 
local coordinates to that of global coordinates; 

Element 1, nodal points Element 2 , nodal points 

Coordinates • Coordinates 

local global local global 



1 1 

2 ‘ 3 

3 4 

4 2 



1 

2 

3 

4 



3 

5 

6 

4 



And noticing that numbering of the generalized displacement, 
q^, follows the relation, 

i = 4m - g (A- 6) 

j = 4n - h (A-7) 

where i and j represent the index for generalized displacements, 

m and n for nodal numbers, g and h are set to be 3 for w, 

2 for w , 1 for w , and 0 for w 

X ' y ' xy 

For example, one of the stiffness matrix elements at 
nodal point 4 in global coordinate can be calculated. 
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b 




a. Nodal points numbering in global coordinate 







b. Nodal points numbering in local coordinate 



Figure 18 Elements and Nodal Points Numbering in 

Global and Local Coordinate for Sample Assembling 




i 
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From equation (A- 6) and (A- 7) , 



i = 16 - g 


(A- 8) 


j = 16 - h 


(A- 9) 



But the same node 4 becomes node 3 at element 1 and node 4 
at element 2 in local coordinate. 

For element 1 , 



(i)j^ = 12 - g 


(A-10) 


(j)j^ = 12 - h 

and for element 2 , 


(A-11) 


,(i )2 = 16 - g 


(A-12) 


(j)2 = 16 - h 


(A-13) 


Therefore, if g = 1 and h = 2, then from equation (A- 
(A-9) , and (A-5) , 


-8) , 


^ij ^ ^15,14 

From equations (A-10) , (A-11) , and (A-1) 


(A- 14) 


11,10 


(A-15) 


and from equations (A-12) , (A- 13) , and (A- 3) , 

^t 

2 15,14 


(A-16) 


Therefore, element of the stiffness matrix in global 
coordinate becomes , 

% * 

^15,14 ^ ^t, '■ '■ ^t^ 

^ 11,10 15,14 


(A-17) 



I 

1 

I 
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All other stiffness matrix element can be expressed by the 
same manner described above, and the programs of assembling 

(13) 

the stxffness matrxx are avaxlable in FEABL. 



I 
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APPENDIX C 

Description of Input Data 

Input data for the programs listed in Appendix D are: 

®o' ^y' ^ ^ described in 

section III. In addition to these, control card for calcula- 
tion of different cases in each job with changing above input 
variables are considered. By using this calculation control 
input card in each job, series of calculation with changing 
some of the input variables can be calculated, and the defor- 
mation behavior with respect to these changing variables can 
easily be visualized by plotting the deformation versus these 
variables . 

Detailed input data card arrangement and format are 



listed in Table 4. 
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Table 4 

Input Data Layout Form 



Card Input Program 



Column Format 


Symbol 


Definition. and Unit 


Card Number 1 


1-5 15 


NCASE 


Number of cases involved 
in each job 


Card Number 2 


1-15 E15.7 


LX 


Length of 1/4 plate in 
X, inches 


16-30 E15.7 


LY 


Length of 1/4 plate in y 


31-35 15 


M 


Number of elements in x 


36-40 15 


N 


Number of elements in y 


Card Niomber 3 


1-12 E12.7 


CX 


Equivalent constrained 
welding moment as in 
Figure 6, Ib-in/in rad. 


13-24 E12.7 


CY 


Equivalent constrained 
welding moment as in 
Figure 6, Ib-in/in rad. 


25-36 E12.7 


THETAZ 


Free joint angular 
change, 9 q, radians 


37-48 E12.7 


GNU 


Poisson's ratio, D 


49-60 E12.7 


THK 


Plate thickness, t, 
inches 


61-72 E12.7 


E 


Young's modulus, E, 

lb/in.2 



Sample 

Input 

14 

16.0E+00 

10. OE+00 
8 
5 

4 . 20E+04 
4 . 20E+04 

55. OE-03 

0.30E+00 
0. 375E+00 

30. OE+06 



APPENDIX D 



Listing of Programs 
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OISTOPTION] A'^ALYSIS OF WFLOf-D PANEL STRUCTURE .... BY D. SHIN 
LX= 0. 16000F 02. LY= 0. lOOOOF 02, M=NO. ElEMt^NTS X PIR= B 
NO. ELEMENTS Y OIR= F, C X= 0.11300F 05, CY= 0.11300E 05 
THFTA 0= 0.45000F-01 RADIANS, POI SSGNS R AT 1 0= 0.30000F 00 
THICKNFSS= 0.37500E 00, YOUNGS MODULUS= 0.30000E 03 
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